{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 17: Signal processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Robert Johansson\n",
    "\n",
    "Source code listings for [Numerical Python - Scientific Computing and Data Science Applications with Numpy, SciPy and Matplotlib](https://www.apress.com/us/book/9781484242452) (ISBN 978-1-484242-45-2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "mpl.rcParams['mathtext.fontset'] = 'stix'\n",
    "mpl.rcParams['font.family'] = 'serif'\n",
    "mpl.rcParams['font.sans-serif'] = 'stix'\n",
    "\n",
    "#sns.set(style=\"whitegrid\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib as mpl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from scipy import fftpack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# this also works:\n",
    "# from numpy import fft as fftpack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from scipy import signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import scipy.io.wavfile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from scipy import io"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spectral analysis of simulated signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def signal_samples(t):\n",
    "    \"\"\" Simulated signal samples \"\"\"\n",
    "    return (2 * np.sin(1 * 2 * np.pi * t) +\n",
    "            3 * np.sin(22 * 2 * np.pi * t) +\n",
    "            2 * np.random.randn(*np.shape(t)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "np.random.seed(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "B = 30.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f_s = 2 * B\n",
    "f_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "delta_f = 0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "N = int(f_s / delta_f)\n",
    "N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "T = N / f_s\n",
    "T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f_s / N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "t = np.linspace(0, T, N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f_t = signal_samples(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(8, 3), sharey=True)\n",
    "axes[0].plot(t, f_t)\n",
    "axes[0].set_xlabel(\"time (s)\")\n",
    "axes[0].set_ylabel(\"signal\")\n",
    "axes[1].plot(t, f_t)\n",
    "axes[1].set_xlim(0, 5)\n",
    "axes[1].set_xlabel(\"time (s)\")\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"ch17-simulated-signal.pdf\")\n",
    "fig.savefig(\"ch17-simulated-signal.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "F = fftpack.fft(f_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f = fftpack.fftfreq(N, 1/f_s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "mask = np.where(f >= 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(3, 1, figsize=(8, 6))\n",
    "axes[0].plot(f[mask], np.log(abs(F[mask])))\n",
    "axes[0].plot(B, 0, 'r*', markersize=10)\n",
    "axes[0].set_xlim(0, 30)\n",
    "axes[0].set_ylabel(\"$\\log(|F|)$\", fontsize=14)\n",
    "\n",
    "axes[1].plot(f[mask], abs(F[mask])/N)\n",
    "axes[1].set_xlim(0, 2)\n",
    "axes[1].set_ylabel(\"$|F|/N$\", fontsize=14)\n",
    "\n",
    "axes[2].plot(f[mask], abs(F[mask])/N)\n",
    "axes[2].set_xlim(19, 23)\n",
    "axes[2].set_xlabel(\"frequency (Hz)\", fontsize=14)\n",
    "axes[2].set_ylabel(\"$|F|/N$\", fontsize=14)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"ch17-simulated-signal-spectrum.pdf\")\n",
    "fig.savefig(\"ch17-simulated-signal-spectrum.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple example of filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "F_filtered = F * (abs(f) <= 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "f_t_filtered = fftpack.ifft(F_filtered)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 3))\n",
    "ax.plot(t, f_t, label='original', alpha=0.5)\n",
    "ax.plot(t, f_t_filtered.real, color=\"red\", lw=3, label='filtered')\n",
    "ax.set_xlim(0, 10)\n",
    "ax.set_xlabel(\"time (s)\")\n",
    "ax.set_ylabel(\"signal\")\n",
    "ax.legend(loc=2)\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"ch17-inverse-fft.pdf\")\n",
    "fig.savefig(\"ch17-inverse-fft.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Windowing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(8, 3))\n",
    "N = 100\n",
    "ax.plot(signal.blackman(N), label=\"Blackman\")\n",
    "ax.plot(signal.hann(N), label=\"Hann\")\n",
    "ax.plot(signal.hamming(N), label=\"Hamming\")\n",
    "ax.plot(signal.gaussian(N, N/5), label=\"Gaussian (std=N/5)\")\n",
    "ax.plot(signal.kaiser(N, 7), label=\"Kaiser (beta=7)\")\n",
    "ax.set_xlabel(\"n\")\n",
    "ax.legend(loc=0)\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"ch17-window-functions.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv('temperature_outdoor_2014.tsv', delimiter=\"\\t\", names=[\"time\", \"temperature\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "df.time = pd.to_datetime(df.time.values, unit=\"s\").tz_localize('UTC').tz_convert('Europe/Stockholm')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df = df.set_index(\"time\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "df = df.resample(\"1H\").ffill()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "df = df[(df.index >= \"2014-01-01\")*(df.index < \"2014-06-01\")].dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df[(df.index >= \"2014-04-01\")*(df.index < \"2014-06-01\")].dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "time = df.index.astype('int')/1e9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "temperature = df.temperature.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "temperature_detrended = signal.detrend(temperature)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "window = signal.blackman(len(temperature_detrended))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temperature_windowed = temperature * window"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "data_fft = fftpack.fft(temperature)\n",
    "data_fft_detrended = fftpack.fft(temperature_detrended)\n",
    "data_fft_windowed = fftpack.fft(temperature_windowed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 4))\n",
    "ax.plot(df.index, temperature, label=\"original\")\n",
    "#ax.plot(df.index, temperature_detrended, label=\"detrended\")\n",
    "ax.plot(df.index, temperature_windowed, label=\"windowed\")\n",
    "ax.set_ylabel(\"temperature\", fontsize=14)\n",
    "ax.legend(loc=0)\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"ch17-temperature-signal.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 4))\n",
    "ax.plot(df.index, temperature_windowed, label=\"original\")\n",
    "ax.plot(df.index, temperature_detrended * window, label=\"windowed\")\n",
    "ax.set_ylabel(\"detrended temperature\", fontsize=14)\n",
    "ax.legend(loc=0)\n",
    "fig.tight_layout()\n",
    "#fig.savefig(\"ch17-temperature-signal.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f = fftpack.fftfreq(len(temperature_windowed), time[1]-time[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "mask = f > 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 4))\n",
    "ax.set_xlim(0.000001, 0.000025)\n",
    "#ax.set_xlim(0.000005, 0.000018)\n",
    "ax.set_xlim(0.000005, 0.00004)\n",
    "\n",
    "ax.axvline(1./86400, color='r', lw=0.5)\n",
    "ax.axvline(2./86400, color='r', lw=0.5)\n",
    "ax.axvline(3./86400, color='r', lw=0.5)\n",
    "ax.plot(f[mask], np.log(abs(data_fft[mask])**2), lw=2, label=\"original\")\n",
    "#ax.plot(f[mask], np.log(abs(data_fft_detrended[mask])**2), lw=2, label=\"detrended\")\n",
    "#ax.plot(f[mask], np.log(abs(data_fft_windowed[mask])**2), lw=2, label=\"windowed\")\n",
    "ax.set_ylabel(\"$\\log|F|$\", fontsize=14)\n",
    "ax.set_xlabel(\"frequency (Hz)\", fontsize=14)\n",
    "ax.legend(loc=0)\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"ch17-temperature-spectrum.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 3))\n",
    "#ax.set_xlim(0.000001, 0.000025)\n",
    "#ax.set_xlim(0.000005, 0.000018)\n",
    "ax.set_xlim(0.000005, 0.00004)\n",
    "\n",
    "ax.axvline(1./86400, color='r', lw=0.5)\n",
    "ax.axvline(2./86400, color='r', lw=0.5)\n",
    "ax.axvline(3./86400, color='r', lw=0.5)\n",
    "\n",
    "y =  np.log(abs(data_fft[mask])**2)\n",
    "ax.plot(f[mask], y / y[10:].max(), lw=1, label=\"original\")\n",
    "\n",
    "y = np.log(abs(data_fft_detrended[mask])**2)\n",
    "ax.plot(f[mask], y / y[10:].max(), lw=2, label=\"detrended\")\n",
    "\n",
    "y = np.log(abs(data_fft_windowed[mask])**2)\n",
    "ax.plot(f[mask], y / y[10:].max(), lw=2, label=\"windowed\")\n",
    "\n",
    "ax.set_ylabel(\"$\\log|F|$\", fontsize=14)\n",
    "ax.set_xlabel(\"frequency (Hz)\", fontsize=14)\n",
    "ax.legend(loc=0)\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"ch17-temperature-spectrum.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spectrogram of Guitar sound"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# https://www.freesound.org/people/guitarguy1985/sounds/52047/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "sample_rate, data = io.wavfile.read(\"guitar.wav\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "sample_rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "data = data.mean(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "data.shape[0] / sample_rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "N = int(sample_rate/2.0); N # half a second"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f = fftpack.fftfreq(N, 1.0/sample_rate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "t = np.linspace(0, 0.5, N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "mask = (f > 0) * (f < 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "subdata = data[:N]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "F = fftpack.fft(subdata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(12, 3))\n",
    "axes[0].plot(t, subdata)\n",
    "axes[0].set_ylabel(\"signal\", fontsize=14)\n",
    "axes[0].set_xlabel(\"time (s)\", fontsize=14)\n",
    "axes[1].plot(f[mask], abs(F[mask]))\n",
    "axes[1].set_ylabel(\"$|F|$\", fontsize=14)\n",
    "axes[1].set_xlabel(\"Frequency (Hz)\", fontsize=14)\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"ch17-guitar-spectrum.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(12, 3))\n",
    "axes[0].plot(t, subdata)\n",
    "axes[0].set_ylabel(\"signal\", fontsize=14)\n",
    "axes[0].set_xlabel(\"time (s)\", fontsize=14)\n",
    "axes[1].plot(f[mask], abs(F[mask]))\n",
    "axes[1].set_ylabel(\"$|F|$\", fontsize=14)\n",
    "axes[1].set_xlabel(\"Frequency (Hz)\", fontsize=14)\n",
    "\n",
    "\n",
    "f_A4 = 440\n",
    "a = 2**(1/12)\n",
    "for note, frequency in [\n",
    "        #('A2', f_A4 * a**(-12-12)),\n",
    "        #('B2', f_A4 * a**(-10-12)),\n",
    "        #('C3', f_A4 * a**(-9-12)),\n",
    "        #('D3', f_A4 * a**(-7-12)),\n",
    "        ('F3', f_A4 * a**(-6-12)),\n",
    "        ('G3', f_A4 * a**(-4-12)),\n",
    "        #('F3', f_A4 * a**(-2-12)),\n",
    "\n",
    "        ('A3', f_A4 * a**(-12)),\n",
    "        #('B3', f_A4 * a**(-10)),\n",
    "        #('C4', f_A4 * a**(-9)),\n",
    "        ('D4', f_A4 * a**(-7)),\n",
    "        #('F4', f_A4 * a**(-6)),\n",
    "        #('G4', f_A4 * a**(-4)),\n",
    "        ('F4', f_A4 * a**(-2)),\n",
    "        ('A4', f_A4),\n",
    "    ]:\n",
    "    axes[1].axvline(frequency, color=\"black\", alpha=0.5)\n",
    "    axes[1].text(frequency*1.01, 2e7, note, fontsize=6)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"ch17-guitar-spectrum.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "N_max = int(data.shape[0] / N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f_values = np.sum(1 * mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "spect_data = np.zeros((N_max, f_values))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "window = signal.blackman(len(subdata))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "for n in range(0, N_max):\n",
    "    subdata = data[(N * n):(N * (n + 1))]\n",
    "    F = fftpack.fft(subdata * window)\n",
    "    spect_data[n, :] = np.log(abs(F[mask]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(8, 6))\n",
    "p = ax.imshow(spect_data, origin='lower',\n",
    "              extent=(0, 1000, 0, data.shape[0] / sample_rate),\n",
    "              aspect='auto',\n",
    "              cmap=mpl.cm.RdBu_r)\n",
    "cb = fig.colorbar(p, ax=ax)\n",
    "cb.set_label(\"$\\log|F|$\", fontsize=16)\n",
    "ax.set_ylabel(\"time (s)\", fontsize=14)\n",
    "ax.set_xlabel(\"Frequency (Hz)\", fontsize=14)\n",
    "\n",
    "\n",
    "\n",
    "f_A4 = 440\n",
    "a = 2**(1/12)\n",
    "for note, frequency in [\n",
    "        #('A2', f_A4 * a**(-12-12)),\n",
    "        #('B2', f_A4 * a**(-10-12)),\n",
    "        #('C3', f_A4 * a**(-9-12)),\n",
    "        ('D3', f_A4 * a**(-7-12)),\n",
    "        #('F3', f_A4 * a**(-6-12)),\n",
    "        ('G3', f_A4 * a**(-4-12)),\n",
    "        ('F3', f_A4 * a**(-2-12)),\n",
    "\n",
    "        ('A3', f_A4 * a**(-12)),\n",
    "        ('B3', f_A4 * a**(-10)),\n",
    "        ('C4', f_A4 * a**(-9)),\n",
    "        ('D4', f_A4 * a**(-7)),\n",
    "        #('F4', f_A4 * a**(-6)),\n",
    "        ('G4', f_A4 * a**(-4)),\n",
    "        ('F4', f_A4 * a**(-2)),\n",
    "        ('A4', f_A4),\n",
    "\n",
    "        ('B4', f_A4 * a**(2)),\n",
    "        ('C5', f_A4 * a**(3)),\n",
    "        ('D5', f_A4 * a**(5)),\n",
    "        ('F5', f_A4 * a**(6)),\n",
    "        ('G5', f_A4 * a**(8)),\n",
    "        ('F5', f_A4 * a**(10)),\n",
    "        ('A5', f_A4 * a**(12)),\n",
    "    ]:\n",
    "    #ax.axvline(frequency, color=\"black\", alpha=0.5)\n",
    "    ax.text(frequency-10, 27, note, fontsize=6)\n",
    "    \n",
    "fig.tight_layout()\n",
    "fig.savefig(\"ch17-spectrogram.pdf\")\n",
    "fig.savefig(\"ch17-spectrogram.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(8, 6))\n",
    "p = ax.imshow(spect_data, origin='lower',\n",
    "              extent=(0, 1000, 0, data.shape[0] / sample_rate),\n",
    "              aspect='auto',\n",
    "              cmap=mpl.cm.RdBu_r)\n",
    "cb = fig.colorbar(p, ax=ax)\n",
    "cb.set_label(\"$\\log|F|$\", fontsize=16)\n",
    "ax.set_ylabel(\"time (s)\", fontsize=14)\n",
    "ax.set_xlabel(\"Frequency (Hz)\", fontsize=14)\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"ch17-spectrogram.pdf\")\n",
    "fig.savefig(\"ch17-spectrogram.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Signal filters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convolution filters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# restore variables from the first example\n",
    "np.random.seed(0)\n",
    "B = 30.0\n",
    "f_s = 2 * B\n",
    "delta_f = 0.01\n",
    "N = int(f_s / delta_f)\n",
    "T = N / f_s\n",
    "t = np.linspace(0, T, N)\n",
    "f_t = signal_samples(t)\n",
    "f = fftpack.fftfreq(N, 1/f_s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "H = (abs(f) < 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "h = fftpack.fftshift(fftpack.ifft(H))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f_t_filtered_conv = signal.convolve(f_t, h, mode='same')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 6))\n",
    "\n",
    "ax = plt.subplot2grid((2,2), (0,0))\n",
    "ax.plot(f, H)\n",
    "ax.set_xlabel(\"frequency (Hz)\")\n",
    "ax.set_ylabel(\"Frequency filter\")\n",
    "ax.set_ylim(0, 1.5)\n",
    "\n",
    "ax = plt.subplot2grid((2,2), (0,1))\n",
    "ax.plot(t - 50, h.real)\n",
    "ax.set_xlabel(\"time (s)\")\n",
    "ax.set_ylabel(\"convolution kernel\")\n",
    "\n",
    "ax = plt.subplot2grid((2,2), (1,0), colspan=2)\n",
    "ax.plot(t, f_t, label='original', alpha=0.25)\n",
    "ax.plot(t, f_t_filtered.real, \"r\", lw=2, label='filtered in frequency domain')\n",
    "ax.plot(t, f_t_filtered_conv.real, 'b--', lw=2, label='filtered with convolution')\n",
    "ax.set_xlim(0, 10)\n",
    "ax.set_xlabel(\"time (s)\")\n",
    "ax.set_ylabel(\"signal\")\n",
    "ax.legend(loc=2)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"ch17-convolution-filter.pdf\")\n",
    "fig.savefig(\"ch17-convolution-filter.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FIR filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "n = 101"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f_s = 1.0 / 3600"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "nyq = f_s/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "b = signal.firwin(n, cutoff=nyq/12, nyq=nyq, window=\"hamming\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "plt.plot(b);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "f, h = signal.freqz(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(8, 3))\n",
    "h_ampl = 20 * np.log10(abs(h))\n",
    "h_phase = np.unwrap(np.angle(h))\n",
    "ax.plot(f/max(f), h_ampl, 'b')\n",
    "ax.set_ylim(-150, 5)\n",
    "ax.set_ylabel('frequency response (dB)', color=\"b\")\n",
    "ax.set_xlabel(r'normalized frequency')\n",
    "ax = ax.twinx()\n",
    "ax.plot(f/max(f), h_phase, 'r')\n",
    "ax.set_ylabel('phase response', color=\"r\")\n",
    "ax.axvline(1.0/12, color=\"black\")\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"ch17-filter-frequency-response.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "temperature_filtered = signal.lfilter(b, 1, temperature)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "temperature_median_filtered = signal.medfilt(temperature, 25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 4))\n",
    "ax.plot(df.index, temperature, label=\"original\", alpha=0.5)\n",
    "ax.plot(df.index, temperature_filtered, color=\"green\", lw=2, label=\"FIR\")\n",
    "ax.plot(df.index, temperature_median_filtered, color=\"red\", lw=2, label=\"median filer\")\n",
    "ax.set_ylabel(\"temperature\", fontsize=14)\n",
    "ax.legend(loc=0)\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"ch17-temperature-signal-fir.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### IIR filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "b, a = signal.butter(2, 14/365.0, btype='high')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "temperature_filtered_iir = signal.lfilter(b, a, temperature)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "temperature_filtered_filtfilt = signal.filtfilt(b, a, temperature)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 3))\n",
    "ax.plot(df.index, temperature, label=\"original\", alpha=0.5)\n",
    "ax.plot(df.index, temperature_filtered_iir, color=\"red\", label=\"IIR filter\")\n",
    "ax.plot(df.index, temperature_filtered_filtfilt, color=\"green\", label=\"filtfilt filtered\")\n",
    "ax.set_ylabel(\"temperature\", fontsize=14)\n",
    "ax.legend(loc=0)\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"ch17-temperature-signal-iir.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# f, h = signal.freqz(b, a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(8, 3))\n",
    "h_ampl = 20 * np.log10(abs(h))\n",
    "h_phase = np.unwrap(np.angle(h))\n",
    "ax.plot(f/max(f)/100, h_ampl, 'b')\n",
    "ax.set_ylabel('frequency response (dB)', color=\"b\")\n",
    "ax.set_xlabel(r'normalized frequency')\n",
    "ax = ax.twinx()\n",
    "ax.plot(f/max(f)/100, h_phase, 'r')\n",
    "ax.set_ylabel('phase response', color=\"r\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filtering Audio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "b = np.zeros(5000)\n",
    "b[0] = b[-1] = 1\n",
    "b /= b.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "data_filt = signal.lfilter(b, 1, data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "io.wavfile.write(\"guitar-echo.wav\",\n",
    "                 sample_rate,\n",
    "                 np.vstack([data_filt, data_filt]).T.astype(np.int16))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# based on: http://nbviewer.ipython.org/gist/Carreau/5507501/the%20sound%20of%20hydrogen.ipynb\n",
    "from IPython.core.display import display\n",
    "from IPython.core.display import HTML\n",
    "def wav_player(filepath):\n",
    "    src = \"\"\"\n",
    "    <audio controls=\"controls\" style=\"width:600px\" >\n",
    "      <source src=\"files/%s\" type=\"audio/wav\" />\n",
    "    </audio>\n",
    "    \"\"\"%(filepath)\n",
    "    display(HTML(src))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "wav_player(\"guitar.wav\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "wav_player(\"guitar-echo.wav\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Versions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%reload_ext version_information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%version_information numpy, matplotlib, scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "npbook_py310",
   "language": "python",
   "name": "npbook_py310"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
